#!/usr/bin/env python2
# -*- coding: utf-8 -*-
"""
#-----------------------------------------------------------------------------
# Name:        tga-noGui.py
# Author:      Cuneyt H. Unlu
# Created:  2013/05/05
# Last mod.:2013/05/21
# Doc:      Conducting thermal analysis using TGA data.
# Prerequisites : Python 2.7, numpy, scipy, matplotlib
#-----------------------------------------------------------------------------
# Changes
# 2013-11-21 Command line interaction was established
#            -p or --pfile argument was ised to process the file
# 2013-12-06 Reaction order (n) was inclined into calculation
"""
import os
import argparse
import re
import math
import time
import numpy
import scipy.signal as sig
import matplotlib.pyplot as plt

class ImportError(Exception):
    """ Import Error class """
    def __init__(self):
        pass
    def __str__(self):
        return repr("File is not imported. Please check field separator")

class CalcTGA:
    """
    Class for estimating activation energy from TGA data.
    """
    def __init__(self):
        """ Initialize class with data supplied from user """
        self.timestamp =  time.strftime("%d-%b-%Y %H:%M:%S",time.localtime())
        self.importData(data["filepath"]+data["filename"],
                        data["field_separator"])
        self.decomposition_begin = data["decomposition_begin"]
        self.decomposition_end = data["decomposition_end"]
        self.calculation_begin = data["calculation_begin"]
        self.calculation_end = data["calculation_end"]
        self.derivative_begin = data["derivative_begin"]
        self.derivative_end = data["derivative_end"]
        self.step = data["step"]
        self.list_breaks = data["list_breaks"]
        self.order = data["order"]

        # convert breaks into centigrade
        self.list_breaks_as_temperature=[]
        for i in self.list_breaks:
            self.list_breaks_as_temperature.append((1000 - 273.15 * i)/i)

        self.findIndex()
        self.calculateDerivative()
        self.prepareSend()
        self.createSend()

    def importData(self,filename,field_separator):
        """ Import data.
            filename : used for import data
            field_separator : to determine columns in csv file
        """
        infile = open(filename, "r")

        self.temperature = []
        self.tg = []
        self.tg_percent = []
        line = infile.readline()

        # check the validity of the field separator
        # if filed separator is wrong there would be
        # an error flag and file will not be imported
        if len(line.split(field_separator)) == 1:
                raise ImportError()
        # instrument filters
        # first line of the ASCII file will be split
        # then first part will be checked
        # regarding to file's origin
        # tg data column, and lines to dismiss will
        # be determined
        if line.split(field_separator)[0] == 'Module':
            self.instrument = 'Pyris'
            dismiss = 20
            tg_column = 3
            factor = 1000.0
        elif line.split(field_separator)[0] == 'CLOSED':
            self.instrument = 'TA Q50'
            dismiss = 42
            tg_column = 2
            factor = 1.0
        else:
            print 'Invalid instrument'

        for line in range(dismiss):
            line = infile.readline()
            if self.instrument == 'Pyris':
                if re.search('Sample Name', line.split(field_separator)[0]):
                    self.sample = str(line.split(field_separator)[1])
            elif self.instrument == 'TA Q50':
                if re.search('Sample', line.split(field_separator)[0]):
                    self.sample = str(line.split(field_separator)[1])
            self.weight(line,factor,field_separator)

        while True:
            line = infile.readline()
            if not line: break
            tmp = line.split(field_separator)
            self.temperature.append(float(tmp[1]))
            self.tg.append(float(tmp[tg_column])/factor)
            self.tg_percent.append(float(tmp[tg_column])/self.W0*100)

        infile.close()

    def weight(self,datum,fac,fs):
        """ find initial weigth """
        tmp = datum.split(fs)
        if self.instrument == 'Pyris':
            if re.search('Sample Weight',tmp[0]):
                self.W0 = float(tmp[1])/fac
        elif self.instrument == 'TA Q50':
            if re.search('Size',tmp[0]):
                self.W0 = float(tmp[1])/fac

    def derivative(self, weight, temperature):
        """ numerical differentiation """
        dif = numpy.diff(weight)/numpy.diff(temperature)
        return dif

    def filt(self,y):
        """ filetr for smoothing derivative plot """
        b,a = sig.butter(4,.01)
        return sig.filtfilt(b,a,y)

    def indices(self,a, func):
        """ Find index values for values in lists """
        return [i for (i, val) in enumerate(a) if func(val)]

    def getIndex (self,a,cond):
        """get index value """
        limitList = self.indices(a, lambda x: x>= cond)
        return limitList[0]

    def stepSize(self,data,size):
        """ find step size as index number """
        i = 1
        while data[i] - data[0] < size:
            i += 1
        return i

    def getData(self,incomingData):
        """ get and assign data to variables """
        data_tg = incomingData[0]
        data_temp = incomingData[1]
        lim1 = incomingData[2]
        lim2 = incomingData[3]
        stp = incomingData[4]
        lim3 = incomingData[5]
        lim4 = incomingData[6]

        x_data = []
        y_data = []
        for i in range(lim3,lim4,stp):
            y = (data_tg [i] - data_tg [lim2]) / (data_tg [lim1] - data_tg [lim2])
            if self.order == 1.0:
                y_element = math.log(1 / y)
            else:
                y_element = (1 - (y**(self.order - 1))) / ((self.order - 1) * (y ** (self.order - 1)))

            x_data.append(1000 / (data_temp [i] + 273.15))
            y_data.append(math.log(y_element))
        return x_data, y_data

    def calculateRegression(self,xdata,ydata):
        """ regression calculation """
        tmp = numpy.polyfit (xdata,ydata,1)
        p = numpy.poly1d(tmp)
        # fit values, and mean
        yhat = [p(z) for z in xdata]
        ybar = sum(ydata)/len(ydata)
        ssreg = sum([ (yihat - ybar)**2 for yihat in yhat])
        sstot = sum([ (yi - ybar)**2 for yi in ydata])
        r = ssreg / sstot
        return -tmp[0] * 8.314, r, yhat

    def findIndex(self):
        """ bind values to variables and get calculation parameters """
        # indices for starting of decomposition
        self.decomposition_begin_index = self.getIndex(
            self.temperature,
            self.decomposition_begin
        )
        self.decomposition_end_index =  self.getIndex(
            self.temperature,
            self.decomposition_end
        )
        # indices for starting of calculation
        self.calculation_begin_index = self.getIndex(
            self.temperature,
            self.calculation_begin
        )
        self.calculation_end_index = self.getIndex(
            self.temperature,
            self.calculation_end
        )
        # indices for derivative calculation
        self.derivative_begin_index = self.getIndex(
            self.temperature,
            self.derivative_begin
        )
        self.derivative_end_index = self.getIndex(
            self.temperature,
            self.derivative_end
        )
        # size of the step as index value
        self.step_index = self.stepSize(
            self.temperature,
            self.step
        )
        # generation of index list for calculation
        self.list_index = [
            self.calculation_begin_index,
            self.calculation_end_index
        ]

        if len(self.list_breaks) > 0:
            for i in self.list_breaks_as_temperature:
                k = self.getIndex(
                    self.temperature,
                    i
                )
                self.list_index.append(k)
        self.list_index.sort()

    def calculateDerivative(self):
        """ calculation of derivative """
        self.range_temperature= self.temperature[
                        self.derivative_begin_index :
                        self.derivative_end_index
            ]
        self.range_tg = self.tg[
                        self.derivative_begin_index :
                        self.derivative_end_index
                    ]
        self.derivative_tg = list(
                            self.derivative(
                                self.range_tg,
                                self.range_temperature
                            )
                        )
        self.derivative_tg_filtered = list(
                            self.filt(
                                self.derivative_tg
                                )
                            )
        # retrieve maximum degradation temperature
        self.maximum_decomposition = self.range_temperature[
                                    self.derivative_tg_filtered.index(
                                        min(self.derivative_tg_filtered)
                                    )
                                ]

    def prepareSend(self):
        """ creates a list to send calculation function """
        self.list_activation_energy = []
        self.list_corr_coeff = []
        self.list_projected_values = []
        self.inverse_temperature = []
        self.inverse_y = []
        self.sent_data = []

        self.sent_data.append(self.tg)
        self.sent_data.append(self.temperature)
        self.sent_data.append(self.decomposition_begin_index)
        self.sent_data.append(self.decomposition_end_index)
        self.sent_data.append(self.step_index)

    def createSend(self):
        """ preparing sending data """
        self.sent_data.append(self.list_index[0])
        self.sent_data.append(self.list_index[1])
        if len(self.list_breaks) == 0:
            output = self.calculateEa(self.sent_data)
            self.getResults(output)
        else:
            for i in range( len( self.list_breaks ) + 1):
                index1 = self.list_index [i]
                index2 = self.list_index[i+1]
                self.sent_data[-2] = index1
                self.sent_data[-1] = index2

                output = self.calculateEa(self.sent_data)

                self.getResults(output)

    def calculateEa(self,coming):
        """ retrieve 1/T, ln ln 1/y as X, Y """
        X, Y = self.getData(coming)
        # calculate activation energy
        activation_energy, corr_coeff, projected_values = self.calculateRegression(X, Y)
        return activation_energy, corr_coeff, projected_values, X, Y

    def getResults(self,output):
        """ gets results from incoming list """
        self.list_activation_energy.append(output[0])
        self.list_corr_coeff.append(output[1])
        self.list_projected_values.append(output[2])
        self.inverse_temperature.append(output[3])
        self.inverse_y.append(output[4])

    def plots(self):
        """ Plots of data. """
        # main figure
        figure = plt.figure(1)
        # inverse_temperature vs. inverse_y,
        plot_inverse = figure.add_subplot(121)
        # tg vs. temperature,
        plot_tga = figure.add_subplot(222)
        # derivative_tg vs. temperature
        plot_derivative = figure.add_subplot(224)

        plot_tga.plot(
            self.temperature,
            self.tg_percent,
            'k-',
            self.temperature[
                self.decomposition_begin_index:
                self.decomposition_end_index
                ],
            self.tg_percent[
                self.decomposition_begin_index:
                self.decomposition_end_index
            ],
            'r-',
            linewidth = 2
        )

        plot_derivative.plot(
            self.temperature[
                self.derivative_begin_index:
                self.derivative_end_index-1
            ],
            self.derivative_tg,
            'k-',
            self.temperature[
                self.derivative_begin_index:
                self.derivative_end_index-1
            ],
            self.derivative_tg_filtered,
            'r-',
            linewidth = 2
        )

        # colors for plot
        # 'k' = black, 'r' = red, 'b' = blue, 'g' = green
        # 'c' = cyan, 'm' = magenta, 'y' = yellow
        colormap = ['k','r','c','m','y','b','g']
        # colored for different sections of plot
        for i in range(len(self.list_breaks)+1):
            plot_inverse.plot(
                self.inverse_temperature[i],
                self.inverse_y[i],
                colormap[i]+'o',
                self.inverse_temperature[i],
                self.list_projected_values[i],
                'y-',
                linewidth = 2
            )

        # Plot formatting
        plot_tga.set_title('Weigth Loss')
        plot_tga.set_xlabel('Temperature (C)')
        plot_tga.set_ylabel('Weigth (%)')
        plot_tga.grid(True)

        plot_derivative.set_title('Derivative Weigth Loss')
        plot_derivative.set_xlabel('Temperature (C)')
        plot_derivative.set_ylabel('Derivative weigth (%/C)')
        plot_derivative.grid(True)

        plot_inverse.set_xlabel('1000/T (1/K)')
        if self.order == 1.0:
            ylabel ='ln ln (1/y)'
        else:
            ylabel ='ln (1-y^(n-1)/(n-1)y^(n-1)'
        plot_inverse.set_ylabel(ylabel)
        plot_inverse.grid(True)

        plt.tight_layout()

        if data["plots"]:
            plt.show()

        if data["report"]:
            plt.savefig(self.sample+'.png', dpi=100, bbox_inches='tight')

    def displayResults(self):
        """ prints results to the console """

        print 40*"="
        print "Date of Analysis:\n\t%s" % self.timestamp
        print "Sample Name:\t%s" % self.sample
        print 40*"="
        print "Limits \t\tBegin\tEnd"
        print 40*"-"
        print "Decomposition:\t%5.2f\t%5.2f" % (
                self.temperature[
                    self.decomposition_begin_index
                    ],
                self.temperature[
                    self.decomposition_end_index
                    ]
                )
        print 40*"-"
        print "Calculation:\t%5.2f\t%5.2f" % (
                self.temperature[
                    self.calculation_begin_index
                    ],
                self.temperature[
                    self.calculation_end_index
                    ]
                )
        print 40*"-"
        print "Derivative:\t%5.2f\t%5.2f" % (
                self.temperature[
                    self.derivative_begin_index
                    ],
                self.temperature[
                    self.derivative_end_index
                    ]

                )
        print 40*"="
        print "Calc. Step Size (C): \t %4.3f" % self.step
        print 40*"="
        print "Max. Degradation (C): \t %5.1f" % self.maximum_decomposition
        print "Reaction Order : \t %2.1f" % self.order
        print "Number of Breaks: \t %d" % len(self.list_breaks)
        print 40*"-"
        if not len(self.list_breaks) == 0:
            print "List of Breaks:\tas C\tas 1/K"
        for i in self.list_index[1:-1]:
            print "\t\t%5.2f\t%6.4f" % (
                    self.temperature[i],
                    1000/(273.15+self.temperature[i]
                        )
                    )
        print 40*"="
        print "Range (C)\tA.E.\tR"
        print 40*"-"
        for i in range(len(self.list_breaks)+1):
            print "%5.2f-%5.2f\t%5.2f\t%5.4f" % (
                    self.temperature[
                        self.list_index[i]
                        ],
                    self.temperature[
                        self.list_index[i+1]
                        ],
                    self.list_activation_energy[i],
                    self.list_corr_coeff[i])
        print 40*"="

    def createReport(self):
        """ Creates an HTML report """

        report_file = open(self.sample+".html",'w')
        report_file.write(
        """
            <!DOCTYPE html>
            <style type="text/css">
            tr:nth-of-type(odd) {
              background-color:#ccc;
            }
            table.table_class {
            	border-width: 0px;
            	border-spacing: 0px;
            	border-style: none;
            	border-collapse: separate;
            }
            table.table_class th {
            	border-width: 0px;
            	padding: 0.2em 2ex 0.2em 2ex;
            	border-style: none;
            }
            table.table_class td {
            	border-width: 0px;
            	padding: 0.2em 2ex 0.2em 2ex;
            	border-style: none;
            	border-color: gray;
            }
            </style>

            <html>
            <head>
            <meta charset="utf-8" />
            <title>Report of the activation energy calculation</title>
            </head>
            <body>
            <h1><b>%s <!-- sample name --></b></h1></br>
            <table>
            <tr>
                <td><h4><em>Time of Analysis</em></h4> </td>
                <td>%s <!-- date and time --></td>
            </tr>
            <tr>
                <td><h4>File path</h4></td>
                <td>%s <!-- file path --></td>
            </tr>
            <tr>
                <td><h4>File name</h4></td>
                <td>%s <!-- file name --></td>
            </tr>
            <tr>
                <td><h4>Instrument </h4></td>
                <td>%s <!-- instrument --></td>
            </tr>
            </table>
            <table class="table_class">
            <tr>
            	<th>Temperature limits (<sup>o</sup>C)</th>
            	<td><b>Begin</b></td>
            	<td><b>End</b></th>
            </tr>
            <tr>
            	<td>Decompostion</td> <!-- Decomposition limits -->
            	<td> %5.2f </td>
                 <td> %5.2f </td>
            </tr>
            <tr>
            	<td>Calculation</td> <!-- Calculation limits -->
            	<td> %5.2f </td>
                 <td> %5.2f </td>
            </tr>
            <tr>
            	<td>Derivative</td> <!-- Derivative limits -->
            	<td> %5.2f </td>
                 <td> %5.2f </td>
            </tr>
            <tr>
            	<td>Calculation step size</td>
            	<td>%4.3f</td>
            	<td></td>
            </tr>
            <tr>
            	<th>Max. degradation</th>
            	<td>%5.2f</td>
            	<td></td>
            </tr>
            <tr>
            	<td><b>Reaction Order</b></td>
            	<td>%2.1f</td>
            	<td></td>
            </tr>
            <tr>
            	<td>Number of breaks</td>
            	<td>%d</td>
            	<td></td>
            </tr>
        """
                %
        (
            self.sample,
            self.timestamp,
            data["filepath"],
            data["filename"],
            self.instrument,
            self.temperature[
                    self.decomposition_begin_index
                    ],
            self.temperature[
                    self.decomposition_end_index
                    ],
            self.temperature[
                     self.calculation_begin_index
                     ],
            self.temperature[
                     self.calculation_end_index
                     ],
            self.temperature[
                     self.derivative_begin_index
                     ],
            self.temperature[
                     self.derivative_end_index
                     ],
            self.step,
            self.maximum_decomposition,
            self.order,
            len(self.list_breaks)
            )

        )

        if not len(self.list_breaks) == 0:
            report_file.write(
                """
                <tr>
                    <td>List of Breaks  </td>
                    <td><sup>o</sup>C   </td>
                    <td>1/K </td>
                </tr>
                """
            )
        for i in self.list_index[1:-1]:
            report_file.write(
                """
                <tr>
                    <td></td>
                    <td> %5.2f </td> <!-- Break point as C -->
                    <td> %6.4f </td> <!-- Break point as 1/K -->
                """
                %
                (
                    self.temperature[i],
                    1000/(273.15 + self.temperature[i])
                )
            )
        report_file.write(
            """
            <tr>
            	<th>Range (C)</th>
            	<th>Act.En.</th>
            	<th>(Corr.Coeff.)</th>
            </tr>
            """
        )

        for i in range(len(self.list_breaks)+1):
            report_file.write(
                """
                <tr>
                    <td>%5.2f &nbsp; - &nbsp; %5.2f </td> <!-- Calculation range -->
                    <td> %5.2f </td> <!-- Activation energy -->
                    <td> %5.4f </td> <!-- Correlation coefficient -->
                </tr>
                """
                %
                    (
                    self.temperature[
                        self.list_index[i]
                        ],
                    self.temperature[
                        self.list_index[i+1]
                        ],
                    self.list_activation_energy[i],
                    self.list_corr_coeff[i]
                    )
                )

        report_file.write(
            """</table>

            <Img src="%s.png" width="600px" alt="" />
            </body>
            </html>
            """
            %
            self.sample
            )

        report_file.close()

    def exportData(self):
        """ exports 1/T and 1/y data """
        export_file = open(self.sample+".XY","w")
        export_file.write("# Inverse temperature and inverse y data\n")
        export_file.write("# Sample %s\n" % self.sample)
        export_file.write("# Date %s \n" % self.timestamp)
        for j in range(len(self.list_breaks)+1):
            for i in range(len(self.inverse_temperature[j-1])):
                export_file.write("%6.5f \t %6.5f\n" % (
                    self.inverse_temperature[j-1][i-1],
                    self.inverse_y[j-1][i-1]
                    )
                )

        export_file.close()

def fileData():
    """ user input part """
    global data
    data = {}
    os.chdir("/home/shared/03-dokuman/00-Belgeler/00-data/TGA/sevim-pva/csv/")
    data["filepath"] = os.getcwd() + os.sep

    parser = argparse.ArgumentParser()

    parser.add_argument('--pfile'   , action='store', dest='pfile')
    parser.add_argument('-p'   , action='store', dest='pfile')

    args = parser.parse_args()

    if args.pfile is None:
        # Argument not given on command line
        args.pfile = "pva.txt"
    execfile(args.pfile)

    # Change data file name
    # sample files
    # PAN homopolymer recorded on TA Q50
#    data["filename"] = "TGA_sample_PAN_Q50.csv"
    # Xylan hemicellulse recorded on TA Q50
    # data["filename"] = "TGA_sample_Xylan_Q50.txt"
    # Xylan hemicellulse recorded on PE Pyris Diamond TG/DTA
    # data["filename"] = "TGA_sample_Xylan_Pyris.txt"

    # Select data file separator, for example \t for tab, etc.
#    data["field_separator"] = ","

    # calculation parameters
    # input section of calculation variables

    # start - end points of decomposition
#    data["decomposition_begin"] = 300
#    data["decomposition_end"] = 350

    # start - end points of calculation
#    data["calculation_begin"] = 305
#    data["calculation_end"] = 320

    # start - end points of derivative calculation
#    data["derivative_begin"] = 300
#    data["derivative_end"] = 350

    # calculation step size in centigrade
#    data["step"] = 0.0001

    # if there are any break point add it to the list
    # value must be 1000/T (Kelvin)
    # values can be directly added from the plot
    # script convert list values to centigrade automatically.
    # i.e.
    # list_breaks = [1.2, 1.4]
    # following notation is also acceptable
    """
    list_breaks = [
        1.2,
        1.4
        ]
    """

#    data["list_breaks"] = [ ]
#    data["order"] = 1.0

    # This section includes Boolean (True/False) values
    # Do you want to see results on screen ?
#    data["screen"] = True
    # Do you want to see plots on screen ?
#    data["plots"] = True
    # Do you want a report as html file ?
#    data["report"] = False
    # Do you want to export data as ASCII files?
#    data["export"] = False

def runCalculation():
    """
    Calculation function of the module.
    """
    global result
    # Don't change anything here
    # if you don't know what you are doing exactly,
    # otherwise script will not work correctly!
    data.update()
    # create Calc instance
    result = CalcTGA()
    result.plots()

    if data["screen"]:
        result.displayResults()

    if data["report"]:
        result.createReport()

    if data["export"]:
        result.exportData()

if __name__ == "__main__":
    # File name and other variables
    fileData()

    # make the calculation
    runCalculation()

# vim: tabstop=8 expandtab shiftwidth=4 softtabstop=4